/*
 * IncompleteLUFactorization.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// An incomplete, level 0, LU factorization preconditioner.
    /// </summary>
    /// <remarks>
    /// The ILU(0) algorithm was taken from: <br/>
    /// Iterative methods for sparse linear systems <br/>
    /// Yousef Saad <br/>
    /// Algorithm is described in Chapter 10, section 10.3.2, page 275 <br/>
    /// </remarks>
    /// <include file='../../../../../examples/trunk/LinearAlgebra/Solvers/Preconditioners/IluFactorization.xml' path='example'/> 
    public sealed class IncompleteLUFactorization : IPreconditioner
    {
        /// <summary>
        /// The original coefficient matrix.
        /// </summary>
        private SparseMatrix mCoefficientMatrix;
        /// <summary>
        /// The matrix holding the lower (L) and upper (U) matrices. The
        /// decomposition matrices are combined to reduce storage.
        /// </summary>
        private SparseMatrix mLuDecomposition;
        /// <summary>
        /// A flag used to indicate if the preconditioner has been created, i.e.
        /// if the <b>mLuDecomposition</b> matrix has been filled with the correct
        /// values.
        /// </summary>
        private bool mWasPreconditionerCreated;

        /// <summary>
        /// Creates an <b>IncompleteLUFactorization</b> instance.
        /// </summary>
        /// <param name="coefficientMatrix">
        /// The <see cref="SparseMatrix"/> upon which the preconditioner is based.
        /// </param>
        /// <exception cref="ArgumentNullException">If <paramref name="coefficientMatrix"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="coefficientMatrix"/> is not a square matrix.</exception>
        public IncompleteLUFactorization(Matrix coefficientMatrix)
        {
            if (coefficientMatrix == null)
            {
                throw new ArgumentNullException("coefficientMatrix", Strings.NullParameterException);
            }

            if (coefficientMatrix.Rows != coefficientMatrix.Columns)
            {
                throw new NotSquareMatrixException();
            }

            mCoefficientMatrix = new SparseMatrix(coefficientMatrix);
        }

        /// <summary>
        /// Returns the decomposed upper triangular matrix.
        /// </summary>
        /// <returns>The upper triangular matrix.</returns>
        internal SparseMatrix UpperTriagonalMatrix()
        {
            return (SparseMatrix)mLuDecomposition.UpperTriangle();
        }

        /// <summary>
        /// Returns the decomposed lower triangular matrix.
        /// </summary>
        /// <returns>The lower triangular matrix.</returns>
        internal SparseMatrix LowerTriagonalMatrix()
        {
            SparseMatrix result = (SparseMatrix)mLuDecomposition.LowerTriangle();
            for (int i = 0; i < result.Rows; i++)
            {
                result[i, i] = 1.0;
            }
            return result;
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        public void CreatePreconditioner()
        {
            mLuDecomposition = new SparseMatrix(mCoefficientMatrix.Rows, mCoefficientMatrix.Columns);

            // M == A
            // for i = 2, ... , n do
            //     for k = 1, .... , i - 1 do
            //         if (i,k) == NZ(Z) then
            //             compute z(i,k) = z(i,k) / z(k,k);
            //             for j = k + 1, ...., n do
            //                 if (i,j) == NZ(Z) then
            //                     compute z(i,j) = z(i,j) - z(i,k) * z(k,j)
            //                 end
            //             end
            //         end
            //     end
            // end
            mCoefficientMatrix.CopyTo(mLuDecomposition);

            for (int i = 0; i < mLuDecomposition.Rows; i++)
            {
                for (int k = 0; k < i; k++)
                {
                    if (mLuDecomposition[i, k] != 0.0)
                    {
                        double t = mLuDecomposition[i, k] / mLuDecomposition[k, k];
                        mLuDecomposition[i, k] = t;
                        if (mLuDecomposition[k, i] != 0.0)
                        {
                            mLuDecomposition[i, i] = mLuDecomposition[i, i] - t * mLuDecomposition[k, i];
                        }
                        for (int j = k + 1; j < mLuDecomposition.Rows; j++)
                        {
                            if (j == i)
                            {
                                continue;
                            }
                            if (mLuDecomposition[i, j] != 0.0)
                            {
                                mLuDecomposition[i, j] = mLuDecomposition[i, j] - t * mLuDecomposition[k, j];
                            }
                        }
                    }
                }
            }

            mWasPreconditionerCreated = true;
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="vector">The right hand side vector.</param>
        /// <param name="result">The left hand side vector.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If <paramref name="vector"/> and <paramref name="result"/> do not
        /// have the same size.</para>
        /// <para> - or -</para>
        /// <para>If the size of <paramref name="result"/> is different the number of rows of the coefficient matrix.</para></exception>
        public void Solve(Vector vector, Vector result)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Count != vector.Count) || (result.Count != mCoefficientMatrix.Rows))
            {
                throw new NotConformableException(Strings.ParametersNotConformable);
            }

            if (!mWasPreconditionerCreated)
            {
                CreatePreconditioner();
            }

            /*
            Solve:
            Lz = y
            Which gives
            for (int i = 1; i < matrix.RowLength; i++)
            {
                z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
            }
            NOTE: l_ii should be 1 because u_ii has to be the value
            */
            Vector rowValues = VectorBuilder.CreateVector(mLuDecomposition.Rows, VectorType.Dense);
            for (int i = 0; i < mLuDecomposition.Rows; i++)
            {
                // Clear the rowValues 
                rowValues.Clear();
                mLuDecomposition.Row(i, rowValues);

                double sum = 0.0;
                for (int j = 0; j < i; j++)
                {
                    sum += rowValues[j] * result[j];
                }

                result[i] = (vector[i] - sum);
            }

            /*    
            Solve:
            Ux = z
            Which gives
            for (int i = matrix.RowLength - 1; i > -1; i--)
            {
                x_i = u_ii^-1 * (z_i - SUM_(j > i) u_ij * x_j)
            }
            */
            for (int i = mLuDecomposition.Rows - 1; i > -1; i--)
            {
                mLuDecomposition.Row(i, rowValues);

                double sum = 0.0;
                for (int j = mLuDecomposition.Rows - 1; j > i; j--)
                {
                    sum += rowValues[j] * result[j];
                }

                result[i] = 1 / rowValues[i] * (result[i] - sum);
            }
        }

        /// <summary>
        /// Solves the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <remarks>
        /// When using the <b>Solve</b> method consider carefully if the
        /// left hand side vector can be reused. If so it would be better to 
        /// use the overloaded <b>Solve</b> method.
        /// </remarks>
        /// <param name="vector">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="vector"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the size of <paramref name="vector"/> is different the number of rows of the coefficient matrix.</para></exception>
        public Vector Solve(Vector vector)
        {
            if (vector == null)
            {
                throw new ArgumentNullException("vector", Strings.NullParameterException);
            }

            if (vector.Count != mCoefficientMatrix.Columns)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            Vector result = VectorBuilder.CreateVector(vector.Count, VectorType.Dense);
            Solve(vector, result);
            return result;
        }
    }
}
